Local Molecular Dynamics with Coulombic Interactions 
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We propose a local, 0{N) molecular dynamics algorithm for the simulation of charged systems. 
The long ranged Coulomb potential is generated by a propagating electric field that obeys modified 
Maxwell equations. On coupling the electrodynamic equations to an external thermostat we show 
that the algorithm produces an effective Coulomb potential between particles. On annealing the 
electrodynamic degrees of freedom the field configuration converges to a solution of the Poisson 
equation much like the electronic degrees of freedom approach the ground state in ab-initio molecular 
dynamics. 



Coulomb's law for the interaction between two charged 
particles is generally presented as a static limit of 
Maxwell's equations valid after all transients have de- 
cayed to zero . Due to the difference between the signal 
propagation speed, c, of electromagnetic radiation and 
excitations in condensed matter, almost all particle-based 
numerical simulations of materials employ the approx- 
imation of static, instantaneous interactions (c = oo). 
This approach has some disadvantages. Since the electro- 
static potential is the unique solution of Poisson's equa- 
tion, even the slightest motion of particles requires a 
global recalculation of the electrostatic potential at ev- 
ery time-step; this calculation can dominate the compu- 
tational effort ^ and represents a major bottleneck for 
the development of efficient multiprocessor codes. 

One might wonder whether more efficient code results 
from a formulation that allows one to reduce the propaga- 
tion speed, but still maintains a sufficiently large separa- 
tion of time scales in a manner familiar from ab-initio 
molecular dynamics 0]. The ratio of the rms parti- 
cle velocity u to c would play the role of an optimiza- 
tion parameter much like the ratio of electron to nu- 
clear masses in quantum chemistry. In order to change c, 
one would simulate the evolution of the coupled particle- 
electroniagnetic system as is routinely done in plasma 
physics y . Such a treatment has the great advantage of 
only requiring local operations, but requires an enormous 
reduction of c in order to be efficient. However, with such 
a dramatic reduction of c the electric field will not fol- 
low the particle motion adiabatically. In this limit there 
is no guarantee that the correct thermodynamic ensem- 
ble is generated. A fundamental question thus arises: Is 
Coulomb's law just a static limit of the Maxwell equa- 
tions or is the law more general? 

Recent work on Monte-Carlo algorithms |^ Q has 
shown that the correct thermodynamic potential is found 
even if the particles and fields propagate at the same 
rate: If one writes that the energy of an electric field E is 



U — j E^/2 d'^r, then a Gibbs distribution characterized 
by interactions in 1/r is generated from the constrained 
integral 

Z({rJ)= /'l?En'5(divE(r)-p({rJ))e-"/'=-^. (1) 

where the charge density, p(r) = '^i^i^i''^ ^ the 
charge of the i'th particle is e^. We work in Heaviside- 
Lorentz electromagnetic units where eg = /io = 1- The 
constraint in the 5- functions is Gauss' law. In the elec- 
trostatic limit, E = — grad 0p, where (j)p is the solu- 
tion of Poisson's equation with source p, but in general 
E — — grad (^p -I- Etr, where Etr is an arbitrary trans- 
verse or rotational vector field. By changing variables 
and integrating over Etr , one immediately sees that 
the longitudinal field components result in a Coulombic 
partition function, while the contribution from the trans- 
verse components merely multiply this partition function 
by a constant. Thus Coulomb's law is valid even in the 
presence of a fully equilibrated and dynamic transverse 
electrostatic field. This result is non-trivial since the re- 
tarded interaction between two charges is not simply 1 Ir 

This Letter implements two molecular dynamics al- 
gorithms which sample Eq. Wc show that our al- 
gorithms generate the correct thermodynamic potential 
and discuss the consequences of lowering c for time de- 
pendent correlations. The algorithms involve only local 
operations on the field degrees of freedom and the CPU 
time per integration step scales linearly with the number 
of particles N. The most obvious choice for sampling 
Eq. 1^ is to directly integrate Maxwell's equations to- 
gether with Newton's equations for the particles coupled 
to the electromagnetic field: 



B = — c curl E, 
E = c curl B - J, 



niiiri = ejE(ri) 

= V, (2) 
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J denotes the electric current due to the particle motion. 
As in the electrostatic limit, we have dropped the Lorentz 
force ei(v X B)/c in the equation for 

Evolution of the system is described by a map in phase 
space. If the Jacobian of this map is unity then there is 
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a conserved density, just like in Hamiltonian mechan- 
ics where LiouviUe's theorem can be apphed: Let us 
integrate Eqs. Q through a time step St and evaluate 
the Jacobian of the transformation dx[/dxj where Xi de- 
notes any one of the variables in Eqs. 13) ■ Since we have 
dx'^/dxi = 1, i.e. the diagonal elements of this matrix 
are unity we find that the Jacobian is 1 -|- 0((5t)^. Thus 
in the limit of small time steps the Jacobian is preserved. 
This implies that a measure conserved by the dynamics 
is 

rf^ = dri^a dVi^a Y[ dEr^p dBr,f3 (3) 
i,a r,/3 

where the products are over the particles and then all 
space. We note that generalized Liouville dynamics 
have turned out very useful recently in the construc- 
tion of new thermostating methods for particle simu- 
lation 0. Since the equations of motion conserve the 
energy U„, = J2 mvf /2 + J d^r + E,^/2} we de- 

duce that the partition function for this system is Z = 
J dfjL 5{Urn — Uq) X (5 (constraints) where the (5-function in- 
cludes all the constraints and conservation law inherent 
in the two Maxwell equations of Eqs. 

A standard hypothesis of ergodicity would lead us to 
guess that we now sample Eq. (^. However, the full 
Maxwell equations have associated with them many in- 
dependent conservation laws in addition to the Gauss 
condition |lO|. most importantly div B = 0. If we simply 
integrate the Maxwell equations we have to supplement 
the Gauss constraint in Eq. |^ with many other con- 
straints in such a way that the analytic formulation of 
the partition function becomes intractable. We should 
reduce the symmetries and conservation laws inherent in 
Maxwell's equations, leaving just Gauss' law; we do this 
by transforming to a constant temperature ensemble and 
coupling the electromagnetic field to thermostats to im- 
prove the ergodicity of the field degrees of freedom. 

We modify two of the equations of motion to 

B = -c curl £-726-1-6, (4) 

where the damping and the noise are related by 
the fluctuation dissipation theorem. The equation for 
the particle velocity is entirely conventional in molecu- 
lar dynamics, that for the magnetic field less so. The 
noise ^ on the magnetic field degrees of freedom is com- 
pletely general; it does not satisfy div ^ ~ 0. Due to 
the coupling of B to the random noise it is ergodic, as is 
Vj. Introduction of the noise has destroyed the unwanted 
constraints arising from Maxwell's equations. However 
taking the divergence of E = c curl B — J we see that 
Gauss' law is still valid for the thermalized equations of 
motion if the equation of continuity div 3 + p ^ holds, 
and we start the simulation with an initial condition con- 
sistent with Gauss' law. 



One can check that the distribution Pq = e~" is a 
fixed point of the thermalized equations. Such a demon- 
stration is often taken as being a sufficient criterion for 
a thermostat in physics; due to the presence of the con- 
servation laws in Maxwell's equations we have examined 
convergence of the dynamics by studying the function 
H = J {P{t) — Pq)'^ / Pq dfi and calculating the dynamics 
of H with a Fokker-Planck equation for the full distribu- 
tion function P. By requiring that H = Q we find that P 
must converge to a function of the form P = A(E, ri)Po 
where 

J d^r (c curl B - J) . A + ^ V,: . ^ (5) 

The condition that this equation is valid for arbitrary 
B leads to the conclusion that ^ is a functional of only 
div E. The general solution of Eq. ©, is then a general 
functional of div E — p: A = A[div E(r) — p({r, Tj})]. 
Choosing div E — p = as a conserved initial condi- 
tion then leads to the required result P{t) Pq. With 
the weight Pq and the measure dfi we integrate over the 
Gaussian variables B and and reproduce the required 
partition function Eq.JQ) independent of c. 

In our implementation of the algorithm, particles of 
mass m move in the continuum. We interpolated the 
charges onto the L^/a^ nodes of a cubic mesh using 3rd 
order B-splines, distributing the charge of each particle 
onto 27 nodes [3 ; higher or lower order schemes are pos- 
sible if they conserve charge. The same mesh is used to 
discretize the field equations, and the electric field E is 
associated with the 'iL^ /a? links. Groups of 4 links are 
grouped into plaquettes, and the magnetic field B lives 
on the 2)L^ /a^ plaquettes The particles interact in 
addition with a shifted Lennard-Jones potential of scale 
a truncated at its minimum, Tc = This discretiza- 

tion is equivalent to a standard 7-point discretization of 
the Laplacian [isj . 

In order to integrate the equations of motion (|2Jl of the 
coupled field/particle system, we use a velocity Verlet 
scheme. First we advance the magnetic field B together 
with the particle velocities to midpoint. Then, the values 
of the B-field and the velocities are used to advance the 
electric field E and the positions r;. Finally, the Langevin 
thermostats Eqs. I@J are applied and the B-field/ velocity 
moves completed. We have used a standard time-step 
5t = O.Olr, where r = \/ma^ jf. is the unit of time and 
t — e^ja the unit of energy. The same time step is used 
for both the particle and field equations. The Courant 
criterion for the stable integration of Maxwell's equations 
is bt < a/^/3c. For the values of c used in this Letter this 
criterion is always satisfied and never limits St. 

The success of the method relies on implementing con- 
straint conserving couplings between particles and fields 
as well as exact local charge continuity. Motion of a par- 
ticle leads to a local (finite) fluctuation of the interpo- 
lated charge density, Api . From this we construct a local 
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FIG. 1: Time for the field integration (□), field-particle cou- 
plings (A) and total time (0) in a system of size L = 30a as 
a function of A'^ on an single AMD Athlon CPU. Also shown 
are results for the same system treated with the P3M method 
using a charge interpolation of the same order (filled symbols, 
■ corresponds to the total reciprocal work). 

current J; such that div J; — —Api. We decompose the 
displacement of a particle Ar^ into a (time reversible) 
sequence of steps {Ax/2, Ay /2, Az, Ay /2, Ax/2} [l3 |. 
Each substep in a direction a leads to a current in only 
those links parallel to a The field update is then 

slaved to the current AE = —A J;. The force acting on 
the particles is found from the principle of virtual work: 
A fluctuation of a particle position da induces a local 
charge fluctuation and thus a local current. The force is 
just fa — —5Um/5a jlQl- This prescription for the force 
leads to the usual electric force /d = ei(E(ri)), where 
(E(ri)) is a local average of the electric field which de- 
pends on the exact form of the interpolation of the charge 
density to the lattice; this is the origin of the acceleration 
in Eq. ©. 

The use of a lattice to discretize the Maxwell equa- 
tions leads to artefacts in the interaction and self en- 
ergy of the particles. These artefacts are removed us- 
ing dynamic subtraction 0. A scalar field that cou- 
ples to the interpolated charges via the energy func- 
tional TyW\ = I [sKVV')^ + mV^] - pV'] d^r leads to 
an effective interaction between particles of the form 
Vy{v) = —eieje~'^'^/r, which when added to the direct 
Coulomb interaction regularizes the short distance sin- 
gularity in l/r. In our molecular dynamics code this 
field obeys the equation of motion 
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(6) 



The Yukawa force fy = +ei{gra.d ip) on the particle 
comes form a local average consistent with the virtual 
work principle, and the total force reads / = /ci + /y + 
/lj. We correct for the Yukawa potential by adding an 
extra analytic Yukawa potential (with opposite sign) to 
the truncated Lennard- Jones potential at short distances. 

Our algorithm runs in two modes. In mode 1, all auxil- 
iary fields are kept at the same temperature as the parti- 
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FIG. 2: Static structure factor S{q) = {s(q)s(— q)), s(q) = 
X^i Si exp (iriq) of a charged symmetric electrolyte {L — 
I5a); three densities n = 0.05a~"^, n — 0.1a~^, n = 0.2a~^ 
and T = e/inkB using the present algorithm in mode I (□), 
c — la/r; transverse electric field and particles thermostated 
to T = e/47rfcs, v/c ^ 0.3 using Eqs. H4I6^ . For ■, the electric 
field is damped, mode II, c = 32a/r, v/c ~ 0.01. Also shown 
is a corresponding Monte-Carlo simulation (A). Solid lines: 
Debye theory. Inset: average instantaneous force for a pair of 
charges, L — 8a for modes I and II. Solid curve: / = —dV/dr, 
where V = —l/Aixr — /QeoL"^ from Ewald summation |I'^ . 
Lattice artefacts removed using dynamic subtraction, ^ = 
0.9/a. 



cles, and we generate the correct thermodynamic interac- 
tion independent of v / c. We also run in "dynamical sim- 
ulated annealing mode" U (mode 11), where we anneal 
the electric degrees of freedom on B to zero temperature 
by removing the noise and using larger values for 72 and 
c, while maintaining the finite temperature thermostat 
on the particles. In this case, our code is closer to tradi- 
tional methods, because the field configuration converges 
to the solution of Poisson's equation when particle mo- 
tion stops. This mode is clearly similar in approach to 
dynamic methods used in quantum simulations _3J with 
v/c a, freely variable optimization parameter. 

Our numerical tests of the algorithm begin with a di- 
rect verification of the 0{N) scaling in Fig. ^ (open sym- 
bols). While the effort to integrate the auxiliary fields 
only depends on the number of grid points, the work re- 
quired for the particle-field couplings rises linearly with 
N . At the "working density" of a typical biomolecular 
simulation u sing a 0.1 nm grid and particle volumes of 
O(0.01nm3) [ill (corresponding to ~ 3000 in Fig.[TJ, 
both parts of the simulation contribute equally to the to- 
tal time. Note in particular that the field integration is 
twice as fast as the reciproca.1 part of a standard Fourier- 
based P3M- implementation [l8j (filled symbols). The ex- 
pense for the particle- field couplings is higher in our code, 
but optimizations in the prefactor can be expected. 

Next, we performed checks on the correctness of the 
method by placing two particles of opposite charges in 
a box and measured the instantaneous force. We com- 
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FIG. 3: Velocity autocorrelation functions for L = 10a, 50 
particles, mode I. Intrinsic damping 71 = 0.1r~^, dotted line: 
exp (— O.lt/r). Solid lines from left to right: correlations for 
c = 0.32,0.35,0.41,0.5,0.7, 1,3.2, lOa/r. Inset: initial decay 
rate ojo from exponential fits. 

pare the results to the short range expansion of the an- 
alytic Ewald summation in the inset of Fig. [3 and find 
excellent agreement for both modes. We also simulated 
a globally neutral electrolyte composed of positive and 
negative particles in order to observe the phenomenon of 
Debye screening. At high enough temperatures, we ex- 
pect a static charge-charge structure factor of the form 
S{q) = e^q^ /{hi^ + q^) where the inverse Debye screen- 
ing length — ne'^/ksT, n is the density of charge 
carriers. Fig. [21 compares S{q) for several densities with 
this expression and also to results measured with our 
related Monte-Carlo algorithm 0. Again we find good 
agreement. As described in we also have studied the 
dispersion law of charge and density fluctuations in a ho- 
mogeneous electrolyte and were able to confirm the fast 
equilibration of density-density and charge-charge corre- 
lations. 

In Fig. [2 we showed that the algorithm generates a 
pair potential equivalent to that of the Ewald summation 
for two very different values of v/c. However, particle dy- 
namics are clearly sensitive to c. As a simple illustration, 
we show velocity-autocorrelations (vi(t)vi(O)) for several 
values of c in Fig.O For very dilute systems the autocor- 
relation is dominated by the damping constant 71 in the 
Langevin thermostat, and (v(i)v(O)) ^ exp(— 71^). At 
higher densities, this decay is modified by collisions, and 
(v(t)v(O)) has a faster decay. The curves in Fig. 13 are 
in this density regime and are sensitive to the intrinsic 
dynamics of the particles. For c < la/r, correlations are 
modified by the low propagation velocities, but saturate 
to a common curve for larger values of c. 

We have simulated a charged system interacting via 
Coulomb forces using an algorithm in which the speed 
of light is a free variable. If v/c is large the dynam- 
ics generate the correct statistical mechanical ensemble 
for Coulomb interacting particles, but dynamic correla- 
tions are modified. On increasing c both the statistical 
mechanical and the local dynamical properties are repro- 



duced correctly. The structure of our code is very simi- 
lar to "particle in cell" plasma codes which are rather 
easy to implement on large multiprocessor computers 
with limited interprocessor bandwidth. We therefore ex- 
pect that on many processors, our algorithm can be com- 
petitive with other fast electrostatic methods including 
Fourier multigrid and fast multipole methods 
[inj i . Our method also generalizes naturally to situations 
with spatial dielectric inhomogeneities, which cannot be 
solved using Fourier techniques and non-standard bound- 
ary conditions, e.g. irregulary shaped volumes. 

We thank Burkhard Diinweg for advice and encourage- 
ment in the implementation of this method. 
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